Water Surface Elevation in 3D¶

In [1]:
import xarray as xr

from schimpy import schism_mesh

import pandas as pd
import holoviews as hv
from holoviews import opts, dim

Open output dataset containing the water surface elevations¶

In [2]:
ds = xr.open_mfdataset('../tests/data/m1_hello_schism/outputs/out2d_*.nc', concat_dim='time', combine="nested",
                  data_vars='minimal', coords='minimal', compat='override')
ds
Out[2]:
<xarray.Dataset>
Dimensions:                  (time: 144, one: 1, nSCHISM_hgrid_node: 2639,
                              nSCHISM_hgrid_face: 4636,
                              nSCHISM_hgrid_edge: 7274,
                              nMaxSCHISM_hgrid_face_nodes: 4, two: 2)
Coordinates:
  * time                     (time) datetime64[ns] 1999-12-31T16:20:00 ... 20...
    SCHISM_hgrid_node_x      (nSCHISM_hgrid_node) float64 dask.array<chunksize=(2639,), meta=np.ndarray>
    SCHISM_hgrid_node_y      (nSCHISM_hgrid_node) float64 dask.array<chunksize=(2639,), meta=np.ndarray>
    SCHISM_hgrid_face_x      (nSCHISM_hgrid_face) float64 dask.array<chunksize=(4636,), meta=np.ndarray>
    SCHISM_hgrid_face_y      (nSCHISM_hgrid_face) float64 dask.array<chunksize=(4636,), meta=np.ndarray>
    SCHISM_hgrid_edge_x      (nSCHISM_hgrid_edge) float64 dask.array<chunksize=(7274,), meta=np.ndarray>
    SCHISM_hgrid_edge_y      (nSCHISM_hgrid_edge) float64 dask.array<chunksize=(7274,), meta=np.ndarray>
Dimensions without coordinates: one, nSCHISM_hgrid_node, nSCHISM_hgrid_face,
                                nSCHISM_hgrid_edge,
                                nMaxSCHISM_hgrid_face_nodes, two
Data variables:
    minimum_depth            (one) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    SCHISM_hgrid             (one) |S1 dask.array<chunksize=(1,), meta=np.ndarray>
    depth                    (nSCHISM_hgrid_node) float32 dask.array<chunksize=(2639,), meta=np.ndarray>
    bottom_index_node        (nSCHISM_hgrid_node) int32 dask.array<chunksize=(2639,), meta=np.ndarray>
    SCHISM_hgrid_face_nodes  (nSCHISM_hgrid_face, nMaxSCHISM_hgrid_face_nodes) float64 dask.array<chunksize=(4636, 4), meta=np.ndarray>
    SCHISM_hgrid_edge_nodes  (nSCHISM_hgrid_edge, two) float64 dask.array<chunksize=(7274, 2), meta=np.ndarray>
    dryFlagNode              (time, nSCHISM_hgrid_node) float32 dask.array<chunksize=(72, 2639), meta=np.ndarray>
    elevation                (time, nSCHISM_hgrid_node) float32 dask.array<chunksize=(72, 2639), meta=np.ndarray>
    depthAverageVelX         (time, nSCHISM_hgrid_node) float32 dask.array<chunksize=(72, 2639), meta=np.ndarray>
    depthAverageVelY         (time, nSCHISM_hgrid_node) float32 dask.array<chunksize=(72, 2639), meta=np.ndarray>
    dryFlagElement           (time, nSCHISM_hgrid_face) float32 dask.array<chunksize=(72, 4636), meta=np.ndarray>
    dryFlagSide              (time, nSCHISM_hgrid_edge) float32 dask.array<chunksize=(72, 7274), meta=np.ndarray>
xarray.Dataset
    • time: 144
    • one: 1
    • nSCHISM_hgrid_node: 2639
    • nSCHISM_hgrid_face: 4636
    • nSCHISM_hgrid_edge: 7274
    • nMaxSCHISM_hgrid_face_nodes: 4
    • two: 2
    • time
      (time)
      datetime64[ns]
      1999-12-31T16:20:00 ... 2000-01-...
      i23d :
      0
      base_date :
      2000 1 1 0.00 8.00
      standard_name :
      time
      axis :
      T
      array(['1999-12-31T16:20:00.000000000', '1999-12-31T16:40:00.000000000',
             '1999-12-31T17:00:00.000000000', '1999-12-31T17:20:00.000000000',
             '1999-12-31T17:40:00.000000000', '1999-12-31T18:00:00.000000000',
             '1999-12-31T18:20:00.000000000', '1999-12-31T18:40:00.000000000',
             '1999-12-31T19:00:00.000000000', '1999-12-31T19:20:00.000000000',
             '1999-12-31T19:40:00.000000000', '1999-12-31T20:00:00.000000000',
             '1999-12-31T20:20:00.000000000', '1999-12-31T20:40:00.000000000',
             '1999-12-31T21:00:00.000000000', '1999-12-31T21:20:00.000000000',
             '1999-12-31T21:40:00.000000000', '1999-12-31T22:00:00.000000000',
             '1999-12-31T22:20:00.000000000', '1999-12-31T22:40:00.000000000',
             '1999-12-31T23:00:00.000000000', '1999-12-31T23:20:00.000000000',
             '1999-12-31T23:40:00.000000000', '2000-01-01T00:00:00.000000000',
             '2000-01-01T00:20:00.000000000', '2000-01-01T00:40:00.000000000',
             '2000-01-01T01:00:00.000000000', '2000-01-01T01:20:00.000000000',
             '2000-01-01T01:40:00.000000000', '2000-01-01T02:00:00.000000000',
             '2000-01-01T02:20:00.000000000', '2000-01-01T02:40:00.000000000',
             '2000-01-01T03:00:00.000000000', '2000-01-01T03:20:00.000000000',
             '2000-01-01T03:40:00.000000000', '2000-01-01T04:00:00.000000000',
             '2000-01-01T04:20:00.000000000', '2000-01-01T04:40:00.000000000',
             '2000-01-01T05:00:00.000000000', '2000-01-01T05:20:00.000000000',
             '2000-01-01T05:40:00.000000000', '2000-01-01T06:00:00.000000000',
             '2000-01-01T06:20:00.000000000', '2000-01-01T06:40:00.000000000',
             '2000-01-01T07:00:00.000000000', '2000-01-01T07:20:00.000000000',
             '2000-01-01T07:40:00.000000000', '2000-01-01T08:00:00.000000000',
             '2000-01-01T08:20:00.000000000', '2000-01-01T08:40:00.000000000',
             '2000-01-01T09:00:00.000000000', '2000-01-01T09:20:00.000000000',
             '2000-01-01T09:40:00.000000000', '2000-01-01T10:00:00.000000000',
             '2000-01-01T10:20:00.000000000', '2000-01-01T10:40:00.000000000',
             '2000-01-01T11:00:00.000000000', '2000-01-01T11:20:00.000000000',
             '2000-01-01T11:40:00.000000000', '2000-01-01T12:00:00.000000000',
             '2000-01-01T12:20:00.000000000', '2000-01-01T12:40:00.000000000',
             '2000-01-01T13:00:00.000000000', '2000-01-01T13:20:00.000000000',
             '2000-01-01T13:40:00.000000000', '2000-01-01T14:00:00.000000000',
             '2000-01-01T14:20:00.000000000', '2000-01-01T14:40:00.000000000',
             '2000-01-01T15:00:00.000000000', '2000-01-01T15:20:00.000000000',
             '2000-01-01T15:40:00.000000000', '2000-01-01T16:00:00.000000000',
             '2000-01-01T16:20:00.000000000', '2000-01-01T16:40:00.000000000',
             '2000-01-01T17:00:00.000000000', '2000-01-01T17:20:00.000000000',
             '2000-01-01T17:40:00.000000000', '2000-01-01T18:00:00.000000000',
             '2000-01-01T18:20:00.000000000', '2000-01-01T18:40:00.000000000',
             '2000-01-01T19:00:00.000000000', '2000-01-01T19:20:00.000000000',
             '2000-01-01T19:40:00.000000000', '2000-01-01T20:00:00.000000000',
             '2000-01-01T20:20:00.000000000', '2000-01-01T20:40:00.000000000',
             '2000-01-01T21:00:00.000000000', '2000-01-01T21:20:00.000000000',
             '2000-01-01T21:40:00.000000000', '2000-01-01T22:00:00.000000000',
             '2000-01-01T22:20:00.000000000', '2000-01-01T22:40:00.000000000',
             '2000-01-01T23:00:00.000000000', '2000-01-01T23:20:00.000000000',
             '2000-01-01T23:40:00.000000000', '2000-01-02T00:00:00.000000000',
             '2000-01-02T00:20:00.000000000', '2000-01-02T00:40:00.000000000',
             '2000-01-02T01:00:00.000000000', '2000-01-02T01:20:00.000000000',
             '2000-01-02T01:40:00.000000000', '2000-01-02T02:00:00.000000000',
             '2000-01-02T02:20:00.000000000', '2000-01-02T02:40:00.000000000',
             '2000-01-02T03:00:00.000000000', '2000-01-02T03:20:00.000000000',
             '2000-01-02T03:40:00.000000000', '2000-01-02T04:00:00.000000000',
             '2000-01-02T04:20:00.000000000', '2000-01-02T04:40:00.000000000',
             '2000-01-02T05:00:00.000000000', '2000-01-02T05:20:00.000000000',
             '2000-01-02T05:40:00.000000000', '2000-01-02T06:00:00.000000000',
             '2000-01-02T06:20:00.000000000', '2000-01-02T06:40:00.000000000',
             '2000-01-02T07:00:00.000000000', '2000-01-02T07:20:00.000000000',
             '2000-01-02T07:40:00.000000000', '2000-01-02T08:00:00.000000000',
             '2000-01-02T08:20:00.000000000', '2000-01-02T08:40:00.000000000',
             '2000-01-02T09:00:00.000000000', '2000-01-02T09:20:00.000000000',
             '2000-01-02T09:40:00.000000000', '2000-01-02T10:00:00.000000000',
             '2000-01-02T10:20:00.000000000', '2000-01-02T10:40:00.000000000',
             '2000-01-02T11:00:00.000000000', '2000-01-02T11:20:00.000000000',
             '2000-01-02T11:40:00.000000000', '2000-01-02T12:00:00.000000000',
             '2000-01-02T12:20:00.000000000', '2000-01-02T12:40:00.000000000',
             '2000-01-02T13:00:00.000000000', '2000-01-02T13:20:00.000000000',
             '2000-01-02T13:40:00.000000000', '2000-01-02T14:00:00.000000000',
             '2000-01-02T14:20:00.000000000', '2000-01-02T14:40:00.000000000',
             '2000-01-02T15:00:00.000000000', '2000-01-02T15:20:00.000000000',
             '2000-01-02T15:40:00.000000000', '2000-01-02T16:00:00.000000000'],
            dtype='datetime64[ns]')
    • SCHISM_hgrid_node_x
      (nSCHISM_hgrid_node)
      float64
      dask.array<chunksize=(2639,), meta=np.ndarray>
      axis :
      X
      location :
      node
      mesh :
      SCHISM_hgrid
      units :
      m
      standard_name :
      projection_x_coordinate
      Array Chunk
      Bytes 20.62 kiB 20.62 kiB
      Shape (2639,) (2639,)
      Dask graph 1 chunks in 2 graph layers
      Data type float64 numpy.ndarray
      2639 1
    • SCHISM_hgrid_node_y
      (nSCHISM_hgrid_node)
      float64
      dask.array<chunksize=(2639,), meta=np.ndarray>
      axis :
      Y
      location :
      node
      mesh :
      SCHISM_hgrid
      units :
      m
      standard_name :
      projection_y_coordinate
      Array Chunk
      Bytes 20.62 kiB 20.62 kiB
      Shape (2639,) (2639,)
      Dask graph 1 chunks in 2 graph layers
      Data type float64 numpy.ndarray
      2639 1
    • SCHISM_hgrid_face_x
      (nSCHISM_hgrid_face)
      float64
      dask.array<chunksize=(4636,), meta=np.ndarray>
      axis :
      X
      location :
      face
      mesh :
      SCHISM_hgrid
      units :
      m
      standard_name :
      projection_x_coordinate
      Array Chunk
      Bytes 36.22 kiB 36.22 kiB
      Shape (4636,) (4636,)
      Dask graph 1 chunks in 2 graph layers
      Data type float64 numpy.ndarray
      4636 1
    • SCHISM_hgrid_face_y
      (nSCHISM_hgrid_face)
      float64
      dask.array<chunksize=(4636,), meta=np.ndarray>
      axis :
      Y
      location :
      face
      mesh :
      SCHISM_hgrid
      units :
      m
      standard_name :
      projection_y_coordinate
      Array Chunk
      Bytes 36.22 kiB 36.22 kiB
      Shape (4636,) (4636,)
      Dask graph 1 chunks in 2 graph layers
      Data type float64 numpy.ndarray
      4636 1
    • SCHISM_hgrid_edge_x
      (nSCHISM_hgrid_edge)
      float64
      dask.array<chunksize=(7274,), meta=np.ndarray>
      axis :
      X
      location :
      edge
      mesh :
      SCHISM_hgrid
      units :
      m
      standard_name :
      projection_x_coordinate
      Array Chunk
      Bytes 56.83 kiB 56.83 kiB
      Shape (7274,) (7274,)
      Dask graph 1 chunks in 2 graph layers
      Data type float64 numpy.ndarray
      7274 1
    • SCHISM_hgrid_edge_y
      (nSCHISM_hgrid_edge)
      float64
      dask.array<chunksize=(7274,), meta=np.ndarray>
      axis :
      Y
      location :
      edge
      mesh :
      SCHISM_hgrid
      units :
      m
      standard_name :
      projection_y_coordinate
      Array Chunk
      Bytes 56.83 kiB 56.83 kiB
      Shape (7274,) (7274,)
      Dask graph 1 chunks in 2 graph layers
      Data type float64 numpy.ndarray
      7274 1
    • minimum_depth
      (one)
      float64
      dask.array<chunksize=(1,), meta=np.ndarray>
      units :
      m
      Array Chunk
      Bytes 8 B 8 B
      Shape (1,) (1,)
      Dask graph 1 chunks in 2 graph layers
      Data type float64 numpy.ndarray
      1 1
    • SCHISM_hgrid
      (one)
      |S1
      dask.array<chunksize=(1,), meta=np.ndarray>
      long_name :
      Topology data of 2d unstructured mesh
      topology_dimension :
      2
      cf_role :
      mesh_topology
      node_coordinates :
      SCHISM_hgrid_node_x SCHISM_hgrid_node_y
      edge_coordinates :
      SCHISM_hgrid_edge_x SCHISM_hgrid_edge_y
      face_coordinates :
      SCHISM_hgrid_face_x SCHISM_hgrid_face_y
      edge_node_connectivity :
      SCHISM_hgrid_edge_nodes
      face_node_connectivity :
      SCHISM_hgrid_face_nodes
      Array Chunk
      Bytes 1 B 1 B
      Shape (1,) (1,)
      Dask graph 1 chunks in 2 graph layers
      Data type |S1 numpy.ndarray
      1 1
    • depth
      (nSCHISM_hgrid_node)
      float32
      dask.array<chunksize=(2639,), meta=np.ndarray>
      units :
      m
      axis :
      Z
      positive :
      down
      location :
      node
      grid_mapping :
      crs
      mesh :
      SCHISM_hgrid
      Array Chunk
      Bytes 10.31 kiB 10.31 kiB
      Shape (2639,) (2639,)
      Dask graph 1 chunks in 2 graph layers
      Data type float32 numpy.ndarray
      2639 1
    • bottom_index_node
      (nSCHISM_hgrid_node)
      int32
      dask.array<chunksize=(2639,), meta=np.ndarray>
      location :
      node
      grid_mapping :
      crs
      mesh :
      SCHISM_hgrid
      Array Chunk
      Bytes 10.31 kiB 10.31 kiB
      Shape (2639,) (2639,)
      Dask graph 1 chunks in 2 graph layers
      Data type int32 numpy.ndarray
      2639 1
    • SCHISM_hgrid_face_nodes
      (nSCHISM_hgrid_face, nMaxSCHISM_hgrid_face_nodes)
      float64
      dask.array<chunksize=(4636, 4), meta=np.ndarray>
      start_index :
      1
      cf_role :
      face_node_connectivity
      Array Chunk
      Bytes 144.88 kiB 144.88 kiB
      Shape (4636, 4) (4636, 4)
      Dask graph 1 chunks in 2 graph layers
      Data type float64 numpy.ndarray
      4 4636
    • SCHISM_hgrid_edge_nodes
      (nSCHISM_hgrid_edge, two)
      float64
      dask.array<chunksize=(7274, 2), meta=np.ndarray>
      start_index :
      1
      cf_role :
      edge_node_connectivity
      Array Chunk
      Bytes 113.66 kiB 113.66 kiB
      Shape (7274, 2) (7274, 2)
      Dask graph 1 chunks in 2 graph layers
      Data type float64 numpy.ndarray
      2 7274
    • dryFlagNode
      (time, nSCHISM_hgrid_node)
      float32
      dask.array<chunksize=(72, 2639), meta=np.ndarray>
      i23d :
      1
      location :
      node
      grid_mapping :
      crs
      mesh :
      SCHISM_hgrid
      Array Chunk
      Bytes 1.45 MiB 742.22 kiB
      Shape (144, 2639) (72, 2639)
      Dask graph 2 chunks in 5 graph layers
      Data type float32 numpy.ndarray
      2639 144
    • elevation
      (time, nSCHISM_hgrid_node)
      float32
      dask.array<chunksize=(72, 2639), meta=np.ndarray>
      i23d :
      1
      location :
      node
      grid_mapping :
      crs
      mesh :
      SCHISM_hgrid
      Array Chunk
      Bytes 1.45 MiB 742.22 kiB
      Shape (144, 2639) (72, 2639)
      Dask graph 2 chunks in 5 graph layers
      Data type float32 numpy.ndarray
      2639 144
    • depthAverageVelX
      (time, nSCHISM_hgrid_node)
      float32
      dask.array<chunksize=(72, 2639), meta=np.ndarray>
      i23d :
      1
      location :
      node
      grid_mapping :
      crs
      mesh :
      SCHISM_hgrid
      Array Chunk
      Bytes 1.45 MiB 742.22 kiB
      Shape (144, 2639) (72, 2639)
      Dask graph 2 chunks in 5 graph layers
      Data type float32 numpy.ndarray
      2639 144
    • depthAverageVelY
      (time, nSCHISM_hgrid_node)
      float32
      dask.array<chunksize=(72, 2639), meta=np.ndarray>
      i23d :
      1
      location :
      node
      grid_mapping :
      crs
      mesh :
      SCHISM_hgrid
      Array Chunk
      Bytes 1.45 MiB 742.22 kiB
      Shape (144, 2639) (72, 2639)
      Dask graph 2 chunks in 5 graph layers
      Data type float32 numpy.ndarray
      2639 144
    • dryFlagElement
      (time, nSCHISM_hgrid_face)
      float32
      dask.array<chunksize=(72, 4636), meta=np.ndarray>
      i23d :
      4
      location :
      face
      grid_mapping :
      crs
      mesh :
      SCHISM_hgrid
      Array Chunk
      Bytes 2.55 MiB 1.27 MiB
      Shape (144, 4636) (72, 4636)
      Dask graph 2 chunks in 5 graph layers
      Data type float32 numpy.ndarray
      4636 144
    • dryFlagSide
      (time, nSCHISM_hgrid_edge)
      float32
      dask.array<chunksize=(72, 7274), meta=np.ndarray>
      i23d :
      7
      location :
      edge
      grid_mapping :
      crs
      mesh :
      SCHISM_hgrid
      Array Chunk
      Bytes 4.00 MiB 2.00 MiB
      Shape (144, 7274) (72, 7274)
      Dask graph 2 chunks in 5 graph layers
      Data type float32 numpy.ndarray
      7274 144
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['1999-12-31 16:20:00', '1999-12-31 16:40:00',
                     '1999-12-31 17:00:00', '1999-12-31 17:20:00',
                     '1999-12-31 17:40:00', '1999-12-31 18:00:00',
                     '1999-12-31 18:20:00', '1999-12-31 18:40:00',
                     '1999-12-31 19:00:00', '1999-12-31 19:20:00',
                     ...
                     '2000-01-02 13:00:00', '2000-01-02 13:20:00',
                     '2000-01-02 13:40:00', '2000-01-02 14:00:00',
                     '2000-01-02 14:20:00', '2000-01-02 14:40:00',
                     '2000-01-02 15:00:00', '2000-01-02 15:20:00',
                     '2000-01-02 15:40:00', '2000-01-02 16:00:00'],
                    dtype='datetime64[ns]', name='time', length=144, freq=None))

Read mesh data of elements and nodes¶

This might not work for quad elements. Will have to deal with those as two non-overlapping triangles

In [3]:
smesh = schism_mesh.read_mesh('../tests/data/m1_hello_schism/hgrid.gr3')

dfelems = pd.DataFrame(smesh.elems,columns=[0,1,2])
#dfelems

dfnodes = pd.DataFrame(smesh.nodes, columns=['x','y','z'])
dfnodes.z = -dfnodes.z
#dfnodes

TriSurface plots from Plotly¶

Needed to add simplices from existing elements information of mesh rather than Delaunay triangulation of nodes

Adapted from https://anaconda.org/philippjfr/brain/notebook?version=2017.05.04.1924¶

The code below allows for the simplices already defined by elems to be used instead of doing a Delaunay triangulation (used from scipy as a way to calculate the simplices)

In [4]:
hv.extension('plotly')

import param

class TriSurface(hv.TriSurface):
    
    simplices = param.Array()

class TriSurfacePlot(hv.plotting.plotly.TriSurfacePlot):

    style_opts = ['cmap', 'plot_edges']

    def get_data(self, element, ranges, style, **kwargs):
        if element.simplices is None:
            return super(TriSurfacePlot, self).get_data(element, ranges, style, **kwargs)
        x, y, z = (element.dimension_values(i) for i in range(3))
        simplices = element.simplices
        return [dict(x=x, y=y, z=z, simplices=simplices)]
    
hv.Store.register({TriSurface: TriSurfacePlot}, 'plotly')
In [5]:
dfelev = dfnodes.copy()

#ds.elevation.min(), ds.elevation.max()

mesh_surface = TriSurface(dfelev, simplices = dfelems.values).opts(plot_edges=False, cmap='gray')
mesh_surface
Out[5]:

Show mesh and water surface¶

The water surface simplices are derived from the nodes with Delaunay triangulation. Can we get the simplices from the information in the schism output files?

In [6]:
dfsurface = dfnodes.copy()

def show_surface(time=0):
    dfsurface.z = ds.elevation.values[time,:]
    water_surface = TriSurface(dfsurface, simplices = dfelems.values).opts(width=800, zlim=(-10,2))
    return water_surface.opts(plot_edges=False, cmap='kbc', clim=(0,2), colorbar=True)
In [7]:
show_surface(0)
Out[7]:

Now combine both to animate both bathymetry and water surface¶

In [8]:
def show_combined(time):
    return show_surface(time)*mesh_surface
In [9]:
import panel as pn
pn.extension()
In [10]:
time_slider = pn.widgets.IntSlider(name='Time Index', start=0, end=len(ds.time)-1)
In [11]:
pn.Row(time_slider, 
       hv.DynamicMap(show_combined, streams={'time':time_slider}).opts(
           width=800, height=800)).servable(title='Water Surface Elevation Animation in 3D')#.show()
Out[11]:
In [ ]: